/***************************************************************************
 *   Copyright (C) 2007 by Alexandre Masselot,,,   *
 *   alex@alex-laptop   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 2 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include "scoringmodeldancik.h"

namespace bol {

  ScoringModelDancik::ScoringModelDancik()
  : SpectraAlignmentScoringModel<PeakTheoMsms, PeakIntensity>()
  {
  }

  ScoringModelDancik::~ScoringModelDancik()
  {
  }

  void ScoringModelDancik::series_add(const string series){
    unsigned int i=PeakTheoMsMsFactory::series2index(series);

    m_series.push_back(i);
  }

  void ScoringModelDancik::init(){
    int n=PeakTheoMsMsFactory::nbSeries();
    m_seriesFlag.assign(n, false);
    for (int i=0; i<m_series.size(); i++){
      // FIXME Alex, please check that - thanks. - Yak -
      // old:  m_seriesFlag[i] =true;
      m_seriesFlag[ m_series[i] ]=true;
    }
    v_score_match.resize(n, 0.L);
    v_score_nomatch.resize(n, 0.L);
    v_bitInfo.resize(n, 0.L);
  }

  double ScoringModelDancik::score(const SpectraAlignment<PeakTheoMsms,PeakIntensity>& spalign){
    double s=0.;
    for(int i=0; i<spalign.index().size(); i++){
      int iseries=(*spalign.spectrumA())[i]->seriesIndex();
      if(spalign.index()[i]>0){
        s+=v_score_match[iseries];
      }else{
        s+=v_score_nomatch[iseries];
      }
    }
    return s;

  }

  void ScoringModelDancik::learn_init(){
    int n=PeakTheoMsMsFactory::nbSeries();
    v_learn_altCpt.assign(n,0);
    v_learn_altNb.assign(n,0);
    v_learn_nullCpt.assign(n,0);
    v_learn_nullNb.assign(n,0);
  }

  void ScoringModelDancik::coverage(const SpectraAlignment<PeakTheoMsms,PeakIntensity>& spalign, unsigned int *nbpot, unsigned int *nbmatch){
    static int lbuf=0;
    static bool **cptPot=NULL;
    static bool **cptMatch=NULL;
    int lmax=0;
    for(vector<PeakTheoMsms*>::const_iterator it=spalign.spectrumA()->begin(); it<spalign.spectrumA()->end(); it++){
      if((*it)->pos() >lmax)
        lmax=(*it)->pos();
    }
    if(lmax>lbuf){
      //we don't really for this reallocation as it will be done 10 times if the largest pept in the bank was to be 40 AA...
      lbuf=lmax+4;
      if(cptPot){
        for(int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
          delete []cptPot[i];
          delete []cptMatch[i];
        }
      }else{
        cptPot=new bool*[PeakTheoMsMsFactory::nbSeries()];
        cptMatch=new bool*[PeakTheoMsMsFactory::nbSeries()];
      }
      for(int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
        cptPot[i]=new bool[lbuf+1];
        cptMatch[i]=new bool[lbuf+1];
      }
    }
    for(int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
      if(m_seriesFlag[i]){
        for (int j=0; j<=lmax; j++){
          cptPot[i][j]=cptMatch[i][j]=false;
        }
      }
    }
    for(int i=0; i<spalign.index().size(); i++){
      unsigned int si=(*spalign.spectrumA())[i]->seriesIndex();
      unsigned int p=(*spalign.spectrumA())[i]->pos();
      if(m_seriesFlag[si]){
        if(spalign.index()[i]>0){
          //there is a match in a regarded series_add
          cptMatch[si][p]=true;
        }
        cptPot[si][p]=true;
      }
    }
    for(int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
      if(m_seriesFlag[i]){
        nbpot[i]=0;
        nbmatch[i]=0;
        for (int j=0; j<=lmax; j++){
          if(cptPot[i][j])
            (*nbpot)++;
          if(cptMatch[i][j])
            (*nbmatch)++;
        }
      }
    }

  }
  /**
   * 
   * @param spalign 
   * @param cpt 
   * @param nb 
   */
  void ScoringModelDancik::learn_push(const SpectraAlignment<PeakTheoMsms,PeakIntensity>& spalign, vector<int>&cpt, vector<int>&nb){
    static unsigned int *nbpot=NULL;
    static unsigned int *nbmatch=NULL;
    if(!nbpot){
      nbpot=new unsigned int[PeakTheoMsMsFactory::nbSeries()];
      nbmatch=new unsigned int[PeakTheoMsMsFactory::nbSeries()];
    }
    coverage(spalign, nbpot, nbmatch);
    for(int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
      if(m_seriesFlag[i]){
        cpt[i]+=nbmatch[i];
        nb[i]+=nbpot[i];
      }
    }
    for(int i=0; i<spalign.index().size(); i++){
      unsigned int si=(*spalign.spectrumA())[i]->seriesIndex();
      if(m_seriesFlag[si]){
        if(spalign.index()[i]>0){
          //there is a match in a regarded series_add
          cpt[si]++;
        }
        nb[si]++;
      }
    }

   }

  void ScoringModelDancik::learn_pushAlternate(const SpectraAlignment<PeakTheoMsms,PeakIntensity>& spalign){
    learn_push(spalign, v_learn_altCpt,v_learn_altNb);
  }
  void ScoringModelDancik::learn_pushNull(const SpectraAlignment<PeakTheoMsms,PeakIntensity>& spalign){
    learn_push(spalign, v_learn_nullCpt,v_learn_nullNb);
  }
  
  void ScoringModelDancik::learn_finalize(){
    int n=PeakTheoMsMsFactory::nbSeries();
    v_score_match.assign(n, 0);
    v_score_nomatch.assign(n, 0);
    v_prob_match.assign(n, 0);
    v_prob_nomatch.assign(n, 0);
    v_bitInfo.assign(n, 0);
    
    for (int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
       if(!m_seriesFlag[i])
         continue;
       cout<<v_learn_altCpt[i]<<"/"<<v_learn_altNb[i]<<"\t"<<v_learn_nullCpt[i]<<"/"<<v_learn_nullNb[i]<<"\n";
       if(v_learn_altNb[i])
         v_prob_match[i]=(1.0*v_learn_altCpt[i])/v_learn_altNb[i];
       if(v_learn_nullNb[i])
         v_prob_nomatch[i]=(1.0*v_learn_nullCpt[i])/v_learn_nullNb[i];
       //we add 10 to avboid the influence of having just one match(or more catastrophic, to have 0)
       double bprob_match=(5.+v_learn_altCpt[i])/(10.+v_learn_altNb[i]);
       double bprob_nomatch=(5.+v_learn_nullCpt[i])/(10.+v_learn_nullNb[i]);
       v_score_match[i]=log(bprob_match)-log(bprob_nomatch);
       v_score_nomatch[i]=log(1.-bprob_match)-log(1.-bprob_nomatch);
       v_bitInfo[i]=bprob_match*v_score_match[i]/log(2.);
    }
  }

  ostream& operator<<(ostream& out, const ScoringModelDancik& model){
    for(int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
      if(!model.m_seriesFlag[i])
        continue;
      out<<PeakTheoMsMsFactory::index2series(i)<<"\tbit="<<model.v_bitInfo[i]
          <<"\tprob_match="<<model.v_prob_match[i]<<"\tprob_nomatch="<<model.v_prob_nomatch[i]
          <<"\tscore_match="<<model.v_score_match[i]<<"\tscore_nomatch="<<model.v_score_nomatch[i]
          <<"\n";
    }
    return out;
  }

  void ScoringModelDancik::setScoreValues(const string series, const double bit, const double score_match, const double score_nomatch){
    for(int i=0; i<PeakTheoMsMsFactory::nbSeries(); i++){
      if(!m_seriesFlag[i])
        continue;
      if(series==PeakTheoMsMsFactory::index2series(i)){
        v_score_match[i]=score_match;
        v_score_nomatch[i]=score_nomatch;
        v_bitInfo[i]=bit;
        return;
      }
    }
    cerr<<"could not ScoringModelDancik::setScoreValues for series "<<series<<endl;
    exit(1);
  }

  void ScoringModelDancik::getScoreValues(const string& series, double* bit, double* score_match, double* score_nomatch) {
    const int idx = PeakTheoMsMsFactory::series2index(series);

    if (m_seriesFlag[idx]) {
      if (bit)		*bit = v_bitInfo[idx];
      if (score_match)	*score_match = v_score_match[idx];
      if (score_nomatch)	*score_nomatch = v_score_nomatch[idx];
    } else {
      if (bit)		*bit = 0.L;
      if (score_match)	*score_match = 0.L;
      if (score_nomatch)	*score_nomatch = 0.L;
    }
  }
}
